A Cellular Automaton Model of Damage 
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We investigate the role of equilibrium methods and stress transfer range in describing the process 
of damage. We find that equilibrium approaches are not applicable to the description of damage 
and the catastrophic failure mechanism if the stress transfer is short ranged. In the long range 
limit, equilibrium methods apply only if the healing mechanism associated with ruptured elements 
is instantaneous. Furthermore we find that the nature of the catastrophic failure depends strongly 
on the stress transfer range. Long range transfer systems have a failure mechanism that resembles 
nucleation. In short range stress transfer systems, the catastrophic failure is a continuous process 
that, in some respects, resembles a critical point. 

PACS numbers: 46.50.+a 62.20.mm 64.60.De 83.60.Uv 



I. INTRODUCTION 

For reasons of both scientific interest and applications 
to materials the subject of damage has interested both 
the physics and material science community for decades. 
Models of damage such as the fiber bundle model (FBM) 
[U, 0] and the hierarchical model Q have been used to ob- 
tain a greater understanding of the mechanisms of dam- 
age and the relationship between damage and phase tran- 
sitions. The treatment of these models has included 
studies of the effect of various ways the individual fibers 
break ^ (t\ , investigations of the nature of failure if heal- 
ing mechanisms are thought to allow the system to re- 
main in equilibrium up to the time of catastrophic failure 
P and the effect of inhomogenieties. [1, Q 

Since these models are often treated with equilibrium 
methods and the time scale of the catastrophic failure 
event makes its description by equilibrium methods in- 
appropriate, we are left to conisder can the approach to 
catastrophic failure be treated with equilibrium meth- 
ods? Equilibrium methods are often applicable, for ex- 
ample, when catastrophic failure is described by nucle- 
ation, where the probability of the occurrence of a critical 
droplet can be calculated by assuming a metastable equi- 
librium state. [m 

In this work we study extensions of the FBM to de- 
termine the applicability of equilibrium methods. In ad- 
dition we probe the nature of the catastrophic failure 
event and how it is affected by the range of stress trans- 
fer. First, we consider the relationship between healing 
and ergodicity which we accomplish using the Thirumalai 
- Mountain (TM) metric, [l^ This will provide informa- 
tion as to the role equilibrium techniques can play in 
descriptions of damage. Second, we look at the nature 
of the critical point that appears to be seen in the global 
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load sharing FBM. 'tI This will naturally bring us to a 
discussion of the role of healing in FBMs. Finally, we 
do a careful analysis of the role the stress transfer range 
has in the nature of the damage process and catastrophic 
failure. 

The structure of the remainder of this paper is as fol- 
lows: In Section [n] we introduce the base model and the 
several variations that we study; in Section [IIII we intro- 
duce the TM metric, describe its application to the FBM 
and our models and discuss the implications of our mea- 
surements on the validity of equilibrium descriptions of 
damage; Q in Section |W] we investigate the nature of 
the critical point in the global load sharing FBM; in 
Section |V] we investigate the impact the stress transfer 
range has on the nature of the damage as well as the 
nature of the catastrophic failure; and in Section IVII we 
present our conclusions. 



II. THE MODEL 

We introduce a continuous, cellular automaton (CA) 
model of damage adapted from the earthquake fault 
model introduced in 1991 by Olami, Feder, and Chris- 
tensen. [Ill The Olami-Feder-Christensen (OFC) model 
is a two-dimensional, CA model motivated by the 
Burridge-Knopoff spring-block model of earthquake 
faults. jl5,] The OFC model is equivalent to the CA 
model proposed by Rundle, Jackson, and Brown (RJB) 
except that in the latter, there is a natural definition 
of internal energy, which makes it simpler to identify 
whether or not the system is in equilibruim. (l6l.ll7| The 
evolution of our model is Markovian and described by 
the following rules. Each site on a lattice (which we take 
as a d = 2 square lattice) is assigned a failure threshold 
CTj- and a residual stress cr[. For the sake of simplicity, in 
this work we will take the failure threshold and residual 
stress to be the same on each site, i.e. a^^^ a^/^ . If a 
site's stress reaches or exceeds its failure threshold, the 
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site reduces its stress to the residual stress by dissipating 
a {ui — cr'') of its stress(where < a < 1 is a parameter 
of the model) and passing the remaining fraction of stress 
(i.e. (1 — ol) {(Ji — a^)) uniformly to the its q ^ neigh- 
bors. The quantity q can range from nearest neighbor 
(g — 4) to "infinite range " where q ^ N is the number 
of sites in the system. We initialize the system by as- 
signing a random stress satisfying a"^ < ai < to each 
site. Given our initializing procedure, it is clear that at 
t = no sites will have ct^ > and hence we must have 
a procedure for inducing failures. We refer to this pro- 
cess as a plate update. There are several ways to do this 
but in this paper we consider the so-called zero velocity 
limit introduced in ref. [3]. According to this proce- 
dure, we search the lattice for the site, i, that minimizes 
— ai. Next, we add an equal amount of stress to each 
site such that the stress on i is now equal to its failure 
threshold. We then discharge the site per the procedure 
above and search the lattice to see if the stress added to 
the neighbors of i caused any of them to fail. If so, we 
discharge their stress as above, and if not, we increase 
the time-step (measured in terms of plate updates) by 
unity and search the lattice for the next site i' , which 
minimizes cr^ — ai. Note that in this version of the model 
a site can still receive and hold stress after it fails. We 
can add noise by resetting sites to a randomized residual 
stress and thus instead of the stress on a sight dropping 
to cr'" it becomes ±rj. This defines the time evolution 
of the OFC model. 

Considerable work has been done on this model with 
noise in the limit that i? — > oo. [H, [l^ll^l The system 
has been shown to be in equilibrium and the probabil- 
ity distribution was shown to be Boltzmann. jl^, [23| 
To better understand the meaning of the results of our 
work, it will be useful to discuss the properties of the 
undamaged model in the R —f oo limit presented in ref- 
erences [ll, [l^ [2^ . We begin with Klein et al [2^1 where 
the authors derived a Langevin equation for the time evo- 
lution of the stress in the RJB model. In this equation, 
all lengths can be scaled by i? oc q^^'''. When this length 
is scaled out, the noise, assumed to by random Gaussian, 
must be scaled according to 

^(f,t)-.^iM. (1) 

In the limit i? — *■ cx3, the Langevin equation becomes lin- 
ear in the stress as all higher order terms are suppressed 
by powers o f 1/ R. This is explained in greater detail in 
Klein et al. [21I [25| In the linearized Langevin equation, 
the drift term can be written as the functional derivative 
of a quadratic action, which guarantees that the prob- 
ability of a state a{\x\/R) approaches a Boltzmann dis- 
tribution as i — » 00. Additionally, by converting the 
Langevin equation into a Fokker-Planck equation where 
the time derivative of the distribution function can be 
written in terms of the divergence of a probability cur- 
rent, it was shown that the stationary solution to this dif- 
ferential equation causes the current to vanish: a general 



definition of equilibrium. By calculating the spectrum of 
eigenvalues of the Fokker-Planck operator, it is scene that 
the Boltzmann solution is the unique, stable solution to 
which all initial conditions evolve in the mean- field limit. 
This was numerically confirmed by Rundle et al 18], in 
which the authors make measurements of the energy in 
the system and find its histogram is given by 

P{E)^g{E)e-f'^, (2) 

where E is the energy stored in the springs, g{E) is the 
density of states which is independent of (3, and f3 is the 
inverse temperature which is related to the amplitude of 
the noise by a fluctuation-dissipation relation. 

As the system is in equilibrium, it is expected that 
the dynamics obey some type of detailed balance. We 
stress that the OFC and RJB models are not in equilib- 
rium for a finite stress transfer range, and so we do not 
expect detailed balance in the non-mean-field case. We 
also note that the Langevin/Fokker-Planck treatment is 
based on coarse graining, which implies that there are 
coarse grained length and time scales below which the 
theory does not describe the system. As such, the dy- 
namics of the model described in the beginning of this 
section do not obey detailed balance as they apply to the 
cells of the automata which exist on a length scale much, 
much smaller than the coarse grained length and all of 
the interactions occur in a single time step which is neces- 
sarily much, much smaller than the coarse graining time. 
Thus, the system can be thought of as in equilibrium, 
only if it is examined on length scales larger than the 
coarse graining size and on time scales greater than the 
coarse graining time and it is at this level that the system 
obeys a type of detailed balance. In the coarse grained 
treatment of the RJB and OFC models the authors take 
the coarse graining length to be the stress transfer range 
R. The coarse graining time is set by the time that the 
system takes to reach its steady state distribution (i.e. 
local equilibrium) in the coarse grained volume. This 
coarse graining time time goes to infinity as the coarse 
grained volume diverges (R 00). These details are 
addressed in full in Ferguson et al. Due to these 

coarse grained length scales, if the system is examined 
on a microscopic level, it may not appear time reversal 
invariant and hence a movie of the system run in reverse 
would look strange. However, if observations are only al- 
lowed at the coarse grained level, then we would see the 
stress fluctuating in an apparently random, and time re- 
versal invariant, manner. Similar effects have been seen 
in references 

In the undamaged OFC or RJB model, there is a 
critical stress, or load, that corresponds to a spinodal, 
^2d\ that marks the limit of the metastable state and is 
responsible for the Gutcnburg-Richter scaling of event 
sizes. [20] A spinodal is a critical point and the scal- 
ing of the avalanche sizes, "earthquake" magnitudes, or 
number of sites that fail in a single plate update, is a con- 
sequence of fluctuations about the spinodal. 0, [2^ The 
theoretical description of the OFC model in the limit of 
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infinite stress transfer range \2d\ has the same physics as 
the theoretical description of the TFB model in Selinger 
et al. Namely they are both described by a Langevin or 
Landau-Ginsburg equation with a one component scalar 
order parameter with the same symmetry. In addition 
the critical slowing down associated with the spinodal as 
calculated from the Langevin equation derived in ref . f20'| 
has the same critical exponent as that associated with the 
time to failure in the global stress transfer TFB models 
studied in ref. One then expects that the behavior 
of the two models is the same in this infinite range stress 
transfer limit. In particular in our model there appears to 
be a metastable state which ends in a spinodal consistent 
with the work of Selinger et al. 

However, the class of FBMs treated by Selinger et al 
and Virgilii et al are unrealistic in that the stress transfer 
range is not infinite in real materials. In addition healing 
of the ruptured elements, if it exists, will not be instan- 
taneous. The studies of the OFC model referenced above 
show that the properties of this class of models depends 
on the stress transfer range. Systems with long but not 
infinite range interactions have pseudo-spinodals, [25l.[26l| 
metastable states with nucleation, and are in a state of 
punctuated equilibrium. That is, they appear to be in 
equilibrium for long periods of time until a large event 
("earthquake") forces the system out of equilibrium. Af- 
ter some relaxation time the system returns to the quasi- 
equilibrium state and the process repeats. Models with 
short range stress transfer (nearest neighbor stress trans- 
fer for example) show no evidence of being in equilibrium 
at any time and also show no evidence of a (pseudo) spin- 
odal. [U, [2^ In this paper we use the OFC model and 
variations that we describe herein to study the effect of 
healing rates, noise, and the stress transfer range. 



A. The Base Model 

The simplest version of our modified model is essen- 
tially the same as the OFC earthquake fault model with 
the difference that after a sites fails a given number of 
times (which we call the site's "number of lives"), it is 
considered dead and no longer interacts with the system. 
Note that this implies that when a site fails within the 
interaction range of a dead site, the live sites receive more 
stress in the transfer process than they would if the site 
were alive. In other words, the stress that would have 
been passed to the dead site is not dissipated, rather, it 
is shared equally among the remaining live sites within 
the interaction range. We have also investigated the case 
in which the stress is passed to dead sites and therefore is 
dissipated. We will discuss the latter case in Section IV. 
But unless otherwise specified, stress is not transfered to 
the dead sites. Given these dynamics, if each site has ten 
lives, then the evolution of the system when no site has 
more than nine failures would be identical to the evolu- 
tion according to the OFC model; however, on the tenth 
failure, the site dies and it no longer holds, or receives. 



stress. In order to get rid of transient effects, the sys- 
tem is run without allowing sites to die (hereafter, we 
call this earthquake mode) for 10® plate updates. After 
the system reaches a steady state we begin evolving it as 
a damage model taking note each time a site fails and 
removing it from the system after a specified number of 
failures. 

Various changes can be made to the base model to ac- 
count for different types of materials. The failure thresh- 
olds can be homogeneous i.e. tr/ — a-^ V z to simulate a 
homogenous or pure material or one could let each site 
have its own failure threshold to better mimic impurities 
in a sample or a heterogeneous material. Another way 
to study the effects of impurities or the general behavior 
of heterogeneous materials would be let each site have a 
different number of lives. 



B. The Step-down Model 

Even if no stress is added to a system the load bearing 
sites weaken over time. For example, a system of fibers 
bearing a constant load will eventually fail. This suggests 
that the failure thresholds in our model should, them- 
selves, be dynamic and decrease over time. We mimic 
this behavior by maintaining a fixed total stress on the 
system and reducing the failure threshold on a site each 
time that site fails. This model can be run in several 
different ways. For example, one could, as in the base 
model, simply specify the number of times a site must fail 
before it dies and reduce the failure threshold by a fixed 
amount each time a site fails. Another method would be 
to define a critical failure threshold such that a site is 
dead once its failure threshold drops below this critical 
value. In addition there is considerable freedom in the 
method of lowering the failure threshold. The threshold 
could be reduced by a random amount or a fixed amount. 
Further, it is known that micro-cracks can heal on 
some time scale. In order to capture this phenomena, the 
failure thresholds could be drawn from a random distri- 
bution whose upper bound decreases in time. Thus, when 
a given site fails, it has some probability of increasing its 
failure threshold, and some probability of decreasing it, 
however, by reducing the upper bound of this distribu- 
tion, one guarantees that, on average, failure thresholds 
will be reduced. In this paper we will only consider step 
down models which kill sites when their failure threshold 
drop below some critical value. Results of simulations of 
heterogeneous systems will be reported in a future pub- 
lication. 



C. The TFB Model 

Since one of our goals is to understand how damage 
affects equilibrium states we will also study the thermo- 
dynamic fiber bundle model (TFB) model introduced by 
Selinger et al. This will serve as a baseline for our stud- 
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ies of the other two classes of models defined above. To 
recover the TFB model and the Disordered Thermody- 
namic Fiber Bundle (DTFB) model of Virgilii et al we 
must take i? — + oo to ensure global load sharing, and we 
must not dissipate stress from the system and hence we 
set a = 0. In both the TFB and DTFB models, the sys- 
tem is described by a Boltzman factor [1, [l^ constructed 
from the hamiltonian 

^ / 1 \ 

J 

where the Sj equals unity for intact fibers and zero for 
broken fibers, Dj is the dissociation energy of the j^^ 
fiber, K is the elastic modulus, and e is the strain on the 
fiber, [i, For the TFB model Dj = D Mj whereas 
for the DTFB model it is fiber dependent [E [3. To 
make contact between the (D)TFB and our models, we 
note that all models discussed herein use the results of 
Brenner [23| to restrict our attention to the Hooke's law 
regime where Ci = nti. In order to simulate these models, 
we use a Metropolis algorithm. 



III. MEASUREMENTS OF ERGODICITY 

As discussed above some theoretical treatments of 
damage utilize equilibrium methods [1, [l3| to obtain an- 
alytic results for simple models of fracture. The assump- 
tion that fracture can be described by an equilibrium 
theory is often justified by the work of Brenner [27| who 
performed experiments on iron whiskers and found that 
the whiskers were, individually, in equilibrium up to the 
point of failure. It is important to note that Brenner's 
concept of equilibrium is not necessarily the same as the 
notion of the equilibrium of statistical ensembles. Bren- 
ner finds that the stress-strain curve of each iron whisker 
displays no hysteresis [l^l up to the point of fracture 
and thus concludes the fibers are in equilibrium. How- 
ever, this is not a sufficient condition for ergodicity of 
systems with many, coupled degrees of freedom. Indeed, 
the process of fracture is an inherently irreversible one 
and hence any model which truly captures the under- 
lying physics must be non-ergodic on some time scale. 
The question of interest is the length of that time scale, 
that is, does the system remain ergodic until the point 
of catastrophic failure or is the physics of damage essen- 
tially non-equilibrium. If it is the latter, then equilibrium 
methods cannot be applied in analytic work. 

In order to test the ergodicity of our model, we mea- 
sure the stress-fiuctuation metric, Q.{t). This metric was 
introduced by Thiumalai and Mountain in 1990 12] and 
adapted to study driven, dissipative systems under stress 
by Ferguson et al fl^. The stress- fluctuation metric is 
a measure of the difference between the time averaged 
stress on a site, CTj(t), and the spacial average of the 
time averaged stress, (^{t)), which approaches the en- 



semble average for ^ 1. Thus, Vl{t) is given by 

m^j^Y^{oAt)-m))\ (4) 

where the sum runs overs the N' non-failed sites on the 
lattice, and the quantities '^j{t) and (o'(t)) are given by 

a,{t) = jj^dt'aj{t'), (5) 

and 

j 

For effectively ergodic systems, n{t) ^ 1/t and hence 
plots of l/n{t) versus t will be linear with positive slope. 

in 

As mentioned above the OFC model has been exhaus- 
tively studied as an earthquake fault model [l^ and we 
know that the model is ergodic provided the interaction is 
long range and some noise is introduced into the system. 
Typically, the noise is added by redrawing the residual 
stress values from a flat, random distribution of width 
Act'" and mean each time a site fails. Thus, if we let 
the long range system evolve for some time in earthquake 
mode before we let the sites die, we know the system will 
be ergodic before the first site dies. Therefore we will be 
able to measure how long the system remains ergodic by 
studying the time-evolution of the metric. 

The first model we test is the TFB model of Selinger et 
al. ]^ We ran the simulation for N = 256^ = 65 536 fibers 
and measured the metric with a/N = 10~^, k = 1, D = 
1, and T = 0.5, following Selinger et al. Q We choose a 
fiber at random, switch its state (e.g. broken to intact) 
and if the energy change is AE < we accept the move 
and if it is AE > we accept the move with probability 
exp {—(3AE). Fibers heal immediately in this model so 
that as soon as a fiber heals, that fiber supports its equal 
share of the load (i.e. its elongation and thus stress is the 
same as the stress on all other intact fibers). The inverse 
temperature f3 is treated as a parameter in the problem. 
As one might expect from dynamics that satisfy detailed 
balance, the TFB model is ergodic (see Fig. 1) since the 
inverse metric is a straight line after some transient time. 
In the work of Selinger et al [8,] the free energy of the TFB 
model is shown theoretically to have a metastable and a 
stable minimum. The failure process is then assumed to 
be a nucleation event for moderate applied global stress. 
However, in our simulations of the TFB model with a 
moderate applied stress the metastable state is has an 
infinite lifetime. This does not mean that the infinitely 
long lived state is not meta-stable in the sense that it is a 
relative minimum of the free energy. Mean-field systems 
in fact are known to have infinitely long lived metastable 
states. [1^ (It should be noted that the simulations in 
Selinger et al are not of the TFB model.) Therefore, the 
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FIG. 1: The inverse TM metric and the order parameter, (j) = 
Nunhioiien/N , (shown in the inset) both as a function of time 
for the TFB with iV = 256^ = 65 536 fibers, a/N = 10"^ 
K = 1, D = 1, and T = 0.5. The metric shows that once 
the order parameter reaches its metastable value, the system 
becomes effectively ergodic. 



TFB vifith global load sharing cannot capture the physics 
of catastrophic failure in a fiber bundle model. '33| 

In addition to the inability of the TFB with global 
load sharing to actually generate the catastrophic fail- 
ure mode, it also incorporates the assumptions that the 
fibers heal instantaneously and the system is in equilib- 
rium. The instantaneous healing assumption is not uni- 
versally applicable and the question of whether or not 
systems with individual fibers that have no hysteresis 
are in equilibrium needs to be tested. By definition of 
the model through a Hamiltonian, the states of the sys- 
tem are described by a Boltzmann distribution and the 
system as a whole is necessarily in equilibrium. Thus, 
to test these equilibrium assumptions, we need to con- 
sider a model where the evolution is described by more 
microscopic physics. This brings us to our set of models. 

First, we examine the base model we introduced in 
Section |TT1 We run these systems with large values of 
the dissipation parameter a to slow down the failure so 
that we may record a significant amount of data as cracks 
appear in the model. We also run the model with three 
different stress transfer ranges, R = 1, R — 30, and with 
R such that g = iV, on a square lattice of size TV with 
periodic boundary conditions, and with Nn^cs = 1, o-^ = 
2.0, a = 0.2, and cr'' = 1.25 ± 0.25, which is the noise as 
described in our description of the model. 

In order that the system be in a steady-state prior to 
any damage, we let the system run in earthquake mode 
for 10^ plate updates. The metrics for the three inter- 
action ranges are shown in Fig. [2] The nearest neighbor 
system is not ergodic even before sites are allowed to 
die. As we can see from the figure, when sites begin to 
die the metric shows an even stronger deviation from er- 
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FIG. 2: The inverse TM metric for the nearest-neighbor 
(R=l) (a), i? = 30 (b), and mean field {q = N) (c) base 
model simulated on a d = 2, periodic, square lattice of linear 
size L = 100,512,768 for R = 1,30, "oo", respectively. The 
parameters of the model are — 2.0, = 1.25 ± 0.25, and 
a = 0.2. The solid line is the inverse TM metric and the 
verticle dashed line indicates the time at which the first site 
dies. 



6 




G 3 
2 



2000 4000 6000 8000 

t 

FIG. 3: The inverse TM metric for the base model when sites 
are healed after one plate updates. The parameters are given 
by J? = 20, L = 512, = 2, cr'' = 1 ± 0.2, and a = 0.1. 
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FIG. 4: The inverse TM metric for the R = 10 step down 
model. The model is simulated on a square lattice of lin- 
ear size L = 256 with periodic boundary conditions and with 
af{t = 0) = 50, ar = 0.25 ± 0.25, and a = 0.05. We con- 
sider a site dead when its failure threshold drops below some 
predefined value which we take as o-^ (t) < + 10~^ A where 
A = cr^(t = 0) - {(j'') = 49.75. 



godicity. Systems with longer range interactions show 
a slightly different behavior. As we know the infinite 
range OFC model is in equilibrium .21j and a finite but 
long range interaction exhibits punctuated equilibrium 
(see the region to the right of the dashed vertical line in 
Figs. |2(b)| and 2(c) I. When sites begin to die, however, 
the systems ceases to be in equilibruim (punctuated or 
otherwise) as is seen in the region after the dashed line 
in Figs. 2(b) and 2(c) We also measured the TM metric 
for the base model with healing and long range {R = 20) 
stress transfer. In this simulation we allow dead sites to 
heal after a proscribed number of time steps. Instanta- 
neous healing is simply the OFC model and we know, 
as mentioned above, that this model is ergodic and inca- 
pable of undergoing catastrophic failure. Clearly, as can 
be seen from the base model, without healing (i.e. the 
healing time oo) the system is not ergodic. In Fig. [3] 
we plot the inverse metric for the base model where the 
dead sites heal after one plate update. As can be seen, 
the system is clearly not ergodic. Measurements for sys- 
tems with longer healing times (not shown) are also not 
ergodic. Thus, the small change from instantaneous heal- 
ing (OFC model) to healing after a single plate update 
(SD model) not only results in a system capable of under- 
going catastrophic failure, it also results in a system that 
is not ergodic on any time interval. Finally, we study the 
step-down (SD) model. The TM metric for the SD model 
with i? = 10 is plotted in Fig. S] Note that in Fig. H the 
data stops near t w 8 x 10^ because the system undergoes 
catastrophic failure in that time step. 



IV. NATURE OF THE CRITICAL POINT 

As we stated in the introduction, the OFC model with 
long range stress transfer has been shown to have a spin- 
odal critical point, [l^l It is the spinodal critical point 
that is responsible for the Gutenburg-Richter scaling in 
the model. Generally in the neighborhood of the spinodal 
the nucleation process is not classical. That is, nucleation 
is not initiated by a compact droplet that has the struc- 
ture of the stable phase in it's interior. Instead, near the 
spinodal, the droplet that initiates nucleation is ramified 
and can be described as a percolation cluster. [2^ [s^l 
As we will see in Section |V] the catastrophic failure mode 
in the long range stress transfer base model appears to 
resemble classical nucleation. This raises the question as 
to how dying sites affect the spinodal seen in the OFC 
model. To answer that question, we let (f> = iVaUvc/A^ 
parameterize the damage in the system, and we run our 
model until = 0.9 (10% of the sites die) and then run 
it in the earthquake mode where 10% of the lattice still 
consists of dead sites, however, we do not kill any addi- 
tional sites. We measure the number of clusters (n^) of 
size s, where a cluster is defined as a set of lattice sites 
that fail as the result of a common "parent" site having 
failed. For example, suppose we force a failure per the 
protocol described in Section HIl and when the forced site 
(the so-called "parent site" ) passes its stress to its neigh- 
bors, three of the neighbors fail. The failed neighbors 
will, in turn, pass their stress to their neighbors. Let 
us further suppose that as a result of one of the three 



7 



neighbors that failed, an additional site fails. Finally, we 
suppose when this site passes its stress, no more sites fail 
and thus the event has stopped (i.e. all sites in the lat- 
tice have Gi < af). In this case, all the sites failed as a 
result of the initial site being forced to fail. These sites, 
including the "parent" site, define a cluster, and in this 
example, the size of the cluster is five. We then run the 
system in the base model mode until (j) = 0.8 and repeat 
the measurement of Ug (not shown). We continue this 
process by similarly decreasing (p. Here we investigate 
two cases: first we consider the case where the dead sites 
still receive stress and thus dissipate it from the system 
and then we consider the case when the dead sites do not 
receive stress at all. 

In Fig. [5] we plot the cluster data on a log-log plot for 
the case in which stress is passed to the dead sites and 
thus dissipated from the system. In the OFC model with 
long range stress transfer, Ug ~ s~^/^ which is consis- 
tent with Fig. 5(a) We can see from Fig. |5(b)] and 5(c) 
that the scaling range decreases and eventually disap- 
pears as (f) decreases from unity. The question remains 
as to whether the motion away from the pseudo-spinodal 
is due to the increased dissipation associated with the 
dead sites or simply due to the dead sites themselves. We 
consider this question below. However, there are two in- 
teresting points associated with the model as run above. 

First, if we consider all of the data generated by the 
various values of (j) and, again, plot the number of clus- 
ters {us) of size s, we get what appears to be a scaling 
law as scene in Fig. [SI The fact that this would appear 
to be a scaling plot is due to the fact that the slope at 
the large cluster end is dominated by the values of (j) near 
one and the contributions from (/> 0.5 are concentrated 
in the region Ug ^ 10. Additionally, we find that run- 
ning the system with these "frozen in" dead sites, seems 
to be the same as running the undamaged model, but 
with a higher dissipation parameter. In fact, the lattice 
described by the damage parameter can be associated 
with an undamaged system running with a dissipation 



a' = 1 — 4>{l — a) , 



(7) 



where a is the dissipation parameter for the system being 
run on a damaged lattice. This is numerically confirmed 
by noting that the two scaling plots generated by the two 
different systems with a and (j) related as above lie one 
on top of the other (see Fig. [7]). 

The second case we consider is when only live neighbors 
and not all of the the neighbors of a failed site evenly 
share the discharged stress. This is the model considered 
in Section [hi] where we plot the TM metric for the OFC 
model with damage. In Fig. [8l we plot the number of 
clusters of size Ug versus s when stress is only transferred 
to live sites. As can be seen from the figure, when the 
stress is transferred only to live sites the scaling is the 
same as the scaling in the pure model. The addition 
of the damage seems to be equivalent to simulating a 
smaller system. This can be seen in Fig. [Sj where we 
compare the scaling plots generated by two systems, one 
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X = 1.49 
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FIG. 5: The number of clusters (ris) of size s plotted on a 
log- log plot for various values of <^ = Ns,uvc/N where stress 
is passed to all site, alive or dead. For (j> ~ 1, we reproduce 
the known results of the OFC model where Us ~ with 
r = 3/2. As we decrease </>, r begins to increase from 3/2 
and the scaling region gets smaller and smaller suggesting 
the damage present in the system drives it away from the 
(pseudo) spinodal. The crosses indicate the data while the 
dashed line is the best fit to the data which gives s ^ 
with r 2.07. 
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FIG. 6: The number of clusters (rig) of size s plotted on a log- 
log plot for all values of <j)- We find that there does appear 
to be a scaling law but the exponent is now approximately 2. 
The crosses indicate the data while the dashed line is the best 
fit to the scaling regime. 



C/3 




S 



FIG. 8: The number of clusters (ua) of size s plotted on a 
log-log plot for various values of <j) where stress is passed only 
to the live sites. For all values of we get the same scaling 
behavior as the pure OFC model where ris ~ with r = 
3/2. 



= 0.8 , a = 0.05 
= 1.0, a = 0.24 
= 0.5 , a = 0.05 
= 1.0, a = 0.525 




0.8 , L = 256 
1 .0 , L = 205 




FIG. 7: The number of clusters (us) of size s plotted on a log- 
log plot where stress is passed to all sites for = 0.8 and (f> = 
0.5 (crosses and exes, respectively) and their corresponding 
undamaged system with a' = 1 — 0(1 — a) (circles for the 
system corresponding to ((> — 0.8 and boxes for the system 
corresponding to = 0.5). 



FIG. 9: The number of clusters (us) of size s plotted on a 
log-log plot where stress is passed only to the live sites for 
</> = 0.8 (exes) and its corresponding undamaged system with 
L' = (boxes). Both scale as ~ with r — 3/2. 



V. LONG V. SHORT RANGE 



of which has linear dimension L and damage parameter 
</)<!, and the other is a non-damaged (0=1) system 
with linear dimension 

L' = ^L. (8) 



In this section, we restrict our attention to the base 
model and study the geometry of the catastrophic fail- 
ure as a function of the stress transfer range. [13] We 
find that the geometry of catastrophic failure in the long 
range stress transfer case is different then in the short 
range case. In the long range model, when the system 
undergoes catastrophic failure the lattice goes from about 
30% dead to 100% dead in one time step (see Fig. [TUl) . 




(c) (d) (c) (d) 



FIG. 10: (Color online) The catastrophic failure event in the 
long range base model. The event begins in a dense region of 
dead sites (a) which, locally, overwhelm the system (b) and 
grow outward (c) in the shape of the interaction, ultimately 
failing the entire lattice (d) in a single time step. 



FIG. 11: (Color online) Two typical nearest-neighbor base 
model lattices at the time of catastrophic failure (percolating 
cluster). Figures (a) and (b) show all of the clusters at the 
time the spanning cluster appears. Figures (c) and (d) isolate 
the percolating cluster from figures (a) and (b), respectively. 



The process begins in a localized region and appears to 
be similar to a nucleation event, with a droplet whose 
interior is the stable phase. In the case of short range 
stress transfer one has to define catastrophic failure a 
bit more carefully. If we are using this model to simu- 
late a FBM then catastrophic failure is defined as 100% 
dead. This is a gradual process which does not resem- 
ble the process in the long range stress transfer model. 
However, we can also think of this model (as well as the 
model with long range stress transfer) as a chip board or 
material such as a rock. In this case, catastrophic failure 
is defined as a cluster of dead sites that span the sys- 
tem or, in other words, the dead sites form a percolating 
cluster. Note that this is a different use of the term clus- 
ter than in Section IIVI Here, cluster refers to a set of 
dead lattice sites that are connected to one and other as 
in nearest-neighbor random site percolation, namely, two 
nearest neighbor dead sites belong to the same cluster. 
In this case, the fraction of dead sites is between 30% 
and 80% at the time of catastrophic failure (percolation) 
and never reaches the state where 100% of the sites are 
dead (see Fig. [TT]) . Note that the critical percolation den- 
sity for nearest-neighbor random site percolation in a two 
dimensional square lattice is approximately 0.593. [2§| 

The short range model has an interesting failure dy- 
namic. As mentioned above, rather than the catastrophic 
failure occurring suddenly, i.e. in one plate update, the 
short range model percolates and reaches a mechanical 



equilibrium (i.e. Ci < a* for all live sites) before it 
reaches the point of 100% failed. The failure begins with 
small independent clusters of dead sites. As the evolu- 
tion continues, some of the clusters begin to merge to 
form larger clusters. Eventually, enough clusters merge 
so that the system is separated into two pieces, one inac- 
cessible to the other without passing through the span- 
ning cluster. Two typical percolating clusters are shown 
in Fig. [TlJ The top two figures show the entire lattice 
with unique clusters corresponding to unique colors and 
the bottom two figures show only the percolating clus- 
ter. We analyzed the fractal dimension d/ of the span- 
ning cluster and found df ~ 1.85 This is consistent with 
the fractal dimension of two dimensional random perco- 
lation, where df = 1.896 [1^, within the accuracy of our 
measurements . 



VI. CONCLUSION 

We have introduced a new model for the study of dam- 
age based on the OFC model for earthquake faults. The 
primary change is that we allow sites to die after a pro- 
scribed number of failures. We study two cases: one in 
which dead sites dissipate all of the stress that is passed 
to them and one in which dead sites are not allowed to 
receive stress. These dead sites mimic damaged elements 
such as broken fibers in the fiber bundle models or cracks 
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in fault systems or materials such as chip boards. Our 
numerical investigation of this model has produced the 
following results: The system ceases to be ergodic, and 
hence is not in equilibrium, in the sense described in Sec- 
tion II, as soon as sites begin to die. Healing the sites 
on a time scale of one plate update or more does not 
restore the ergodicity seen in the OFC model. The pres- 
ence of dead sites seems to not only drive the system out 
of equilibrium but also drives it away from the (pseudo) 
spinodal in the case of long range stress transfer where 
the stress is transferred to the dead sites causing addi- 
tional dissipation. If the stress is only transferred to live 
sites then the system remains near the (pseudo) spinodal. 

Catastrophic failure in the long range system resem- 
bles classical nucleation. However, this requires further 
investigation. First, the system is not in metastable equi- 
librium as can be seen from the TM metric data so the 
standard quasi-equilibrium methods [Tol [Tl| do not ap- 
ply. Second, the "nucleation" process that causes catas- 
trophic failure is not seen in the nearest neighbor stress 



transfer system so there exists some crossover regime that 
needs to be studied. Catastrophic failure in the short 
range stress transfer system does not resemble nucleation 
but can be classified as a continuous process. In the case 
where catastrophic failure is defined as the lattice being 
split into two separated pieces by a percolating cluster 
of dead sites, analogous to the fracturing of a chip board 
[S^l, the catastrophic failure event can be classified as a 
fractal. The research presented on this model suggests 
several directions for further investigations into the na- 
ture of damage and catastrophic failure. 
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